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Abstract Functionals of Brownian motion have diverse applications in physics, 
mathematics, and other fields. The probability density function (PDF) of 
Brownian functionals satisfies the Feynman-Kac formula, which is a Schrodinger 
equation in imaginary time. In recent years there is a growing interest in par- 
ticular functionals of non-Brownian motion, or anomalous diffusion, but no 
equation existed for their PDF. Here, we derive a fractional generalization 
of the Feynman-Kac equation for functionals of anomalous paths based on 
sub-diffusive continuous-time random walk. We also derive a backward equa- 
tion and a generalization to Levy flights. Solutions are presented for a wide 
number of applications including the occupation time in half space and in an 
interval, the first passage time, the maximal displacement, and the hitting 
probability. We briefly discuss other fractional Schrodinger equations that 
recently appeared in the literature. 



1 Introduction 

A Brownian functional is defined as A = J Q * U[x(t)]cIt, where x(t) is a trajec- 
tory of a Brownian particle and U(x) is a prescribed function [TJ. Functionals 
of diffusive motion arise in numerous problems across a variety of scientific 
fields from condensed matter physics (2,3,4 , to hydrodynamics [5], meteo- 
rology 6:, and finance [TIE]. The distribution of these functionals satisfies 
a Schrddinger-like equation, derived in 1949 by Kac inspired by Feynman's 
path integrals [9 . Denote by G(x, A, t) the joint probability density function 
(PDF) of finding, at time t, the particle at x and the functional at A. The 
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Feynman-Kac theory asserts that (for U(x) > 0) [TJ|9] 

^G(x,p,t) = K^-G{x,p,t) -pU(x)G(x,p,t), (1) 

where the equation is in Laplace space, A — > p, and K is the diffusion coef- 
ficient. 

The celebrated Feynman-Kac equation (JTJ) describes functionals of normal 
Brownian motion. However, we know today that in a vast number of systems 
the underlying processes exhibit anomalous, non-Brownian sub-diffusion, as 
reflected by the nonlinear relation: (x 2 ) ~ t a , < a < 1 [T0llllU121ll3lfl4] . 
While a few specific functionals of anomalous paths have been investigated 
[TS1[TB] , a general theory is still missing. 

Several functionals of anomalous diffusion arc of interest. For example, the 
time spent by a particle in a given domain, or the occupation time, is given 
by the functional A = L U [x(t)]cLt, where U(x) = 1 in the domain and is 
zero otherwise [171,18, 19,2(7]. Such a functional can be used in kinetic studies 
of chemical reactions that take place exclusively in the domain. Consider 
for example a particle diffusing in a medium containing an interval that 
is absorbing at rate R. The average survival probability of the particle is 
(exp(— RA)) [2T]. Two other related functionals are the occupation time in 
the positive half-space (U(x) = 0(x)) and the local time (U(x) = 5(x)) [T5] 

Another interesting family of functionals arises in the study of NMR 
[25) . In a typical NMR experiment, the macroscopic measured signal can be 
written as E = (e lv ) where ip = 7 J„ B{x(t)]cIt is the phase accumulated 
by each spin, 7 is the gyromagnetic ratio, B(x) is a spatially-inhomogeneous 
external magnetic field, and x(t) is the trajectory of each particle. NMR 
therefore indirectly encodes information regarding the motion of the particles. 
Common choices of the magnetic field B are B(x) — x and B(x) = x 2 [25] . 
For dispersive systems with inhomogeneous disorder where the motion of the 
particles is non-Brownian, the phase (p is a non-Brownian functional with 
U(x) — x or U(x) — x 2 . 

In this paper, we develop a general theory of non-Brownian function- 
als. The process we consider as the mechanism that leads to non-Brownian 
transport is the sub-diffusive continuous-time random- walk (CTRW). This is 
an important and widely investigated process that is frequently used to de- 
scribe the motion of particles in disordered systems [TtJllTTl[T2 H 2 6 [|2"T] . In the 
scaling limit of this process, we derive the following fractional Feynman-Kac 
equation: 

^G(x,p,t) = K a -^Vl' a G(x,p,t) -pU(x)G(x,p,t), (2) 

where the symbol T>\~ a is Friedrich's substantial fractional derivative and is 
equal in Laplace space t — >• s to [s+pU(x)] 1 ~ a [35]. In the rest of the paper, 
we derive Eq. ^ and its backward version and then investigate applications 
for specific functionals of interest. A brief report of part of the results has 
recently appeared in [29] . 
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2 Derivation of the equations 

We use the continuous-time random-walk (CTRW) model as the underlying 
process leading to anomalous diffusion [T0llll[ll2[l26ll27] . In CTRW, an infi- 
nite one-dimensional lattice with spacing a is assumed, and allowed jumps 
are to nearest neighbors only and with equal probability of jumping left or 
right. Waiting times between jump events are independent identically dis- 
tributed random variables with PDF ip(r), and the process starts with a 
particle at x = Xq. The particle waits at xq for time r drawn from ip(r) 
and then jumps with probability 1/2 to either xq + a or x$ — a, after which 
the process is renewed. We assume that no external forces are applied and 
that for long waiting times, ip[r) ~ B a r' (1+a) /\r(-a)\. For < a < 1, 
the average waiting time is infinite and the process is sub-diffusive with 
(x 2 ) = 2K a t a /r{\ + a) (K a = a 2 /(2B a ), units m 2 /sec Q ) [3D]. We look for 
the differential equation that describes the distribution of functionals in the 
scaling limit of this model. 



2.1 Derivation of the fractional Feynman-Kac equation 

Recall that the functional is defined as A = f Q U{x(t)]cIt and that G(x, A, t) 
is the joint PDF of x and A at time t. For the particle to be at (x, A) at 
time t, it must have been at [x, A — tU(x)] at the time t — r when the last 
jump was made. Let Q n (x, A, t)dt be the probability of the particle to make 
its nth jump into (x, A) in the time interval [t, t + dt]. Thus, 

r-t OO 

G(x,A,t)= / W(T)y2Q n [x,A-TU(x),t-T]dr, (3) 
J o „ =0 

where W(t) = 1 — J q t ip{ T ')dT' is the probability for not moving in a time 
interval of length r. 

To arrive into (x,A) after n + 1 jumps, the particle must have arrived 
after n jumps into either [x — a, A — tU(x — a)] or [x + a, A — tU (x + a)} , 
where r the time between the jumps. Since the probabilities of jumping left 
and right are equal, we can write a recursion relation for Q n : 

Q n+1 (x,A,t) = J x/}(r) ^~Q n [x + a,A-TU(x + a),t-T] (4) 
+ ^Q n [x — a,A — tU(x — a),t — r]^ dr, 

where i/j(t) is the PDF of r, the time between jumps. For n — (no jumps 
were made), Qo — S(x — xo)S(A)S(t). 

Assume that U (x) > for all x and thus A > (an assumption we will 
relax later). Let Q n (x,p, t) be the Laplace transform A — > p of Q n (x, A, t) — 
we use along this work the convention that the variables in parenthesis define 
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the space we are working in. We note that 



e- pA Q n [x,A- rU{x),t]dA = e - pTU{x) / e~ pA ' Q n (x, A' , t)dA' 



where we used the fact that Q n (x, A,t) = for A < 0. Thus, Laplace trans- 
forming A — > p Eq. Q we find 



i ft 

Q n+1 (x,p,t) = - / ^(T)e~ pTU ^Q n {x + a lP ,t-T)dT 



+ \ f HT)e~ pTU[x ~ a) Qn(x-a, Pl t~T)dT. (5) 



Laplace transforming t — > s Eq. ([5]) using the convolution theorem, 



l+ l(x,p,s) = -ijj[s + pU(x + a)]Q n (x + a,p, s) 
+ '~i>[s+pU(x - a)]Q n (x - a,p,s), 



(6) 



where ip(s) is the Laplace transform of the waiting time PDF. Fourier trans- 
forming x — > k Eq. (jnj), 



n+ i(fc,p, s) = cos(fca) / e ifc:r -!/'[s + pU(x)]Q n (x,p, s)dx. 



Applying the Fourier transform identity J-{xf(x)} = —i-Srf(k) 



Q n+ i(k,p, s) = cos(fca)-0 



Q n (k,p,s) 



(7) 



Note that the order of the terms is important: tp [s + pU (— i-^)] does not 
commute with cos(fca). Summing Eq. (JJJ over all n, using the initial condition 
Qo(k,p, s) = e lkx ° , and rearranging, we obtain, 



Q n (k,p,s) = < 1 — cos(ka)'ip 

« — n v 



S + pC/ ( — i 



P 



(8) 



We next use our expression for Xl^Lo *9™ ^° calculate G(x, A, t). Trans- 
forming Eq. (]3]) (x, A,t) — ► (k,p,s), 



G(k 7 p, S) = 1 ( d , dk>l Qn(k, P , s), 



(9) 
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where we used W(s) = J °° e st 1 — J* * ip(r)dT dt = [1 — ?/>(s)]/s. Substitut- 
ing Eq. ([8]) into Eq. ([9]), we find the formal solution 



G(k,p, s) 



x < 1 — cos(fca)V> 



(10) 



To derive a differential equation for G(x^p,t), we recall the waiting time 
distribution is ip(t) ~ £? Q i~( 1+Q *'/|/ n (— a)| and write its Laplace transform 
V>(s) for s -> as [12] 



1 - B a s a 



< a < 1 , s ->• 0. 



(11) 



Substituting Eq. (fTTj) into Eq. (|f Op . applying the small fc expansion cos(fca) 
f — fc 2 a 2 /2, and neglecting the high order terms, we have 



G{k,p,s) 



s + pU —i 



d_ 

dk 



k 



s + pU —i 



d_ 

dk 



,ikx 



where we used the generalized diffusion coefficient K a = lim a 2 s a ->o a2 /(2S Q ) 
[30| . By neglecting the high order terms in s and k we effectively reach the 
scaling limit of the lattice walk [3"T1I32"II53"] . Rearranging the expression in the 
last equation we find 



sG(k,p, s) 



K a k 2 
pU 



'dk 



s + pU 

G(k,p,s) 



G(k,p, s) 



Inverting k — > x, s — > t we finally obtain our fractional Feynman-Kac equa- 
tion 

i-G(x,p,t) = K a ^V]- a G{x,p,t) -pU(x)G(x,p,t). (12) 



dt 



dx 2 



The initial condition is G(x,A,t = 0) = S(x — Xq)S(A), or G(x,p,t = 0) 
S(x - x Q ). 



T> t a is the fractional substantial derivative operator introduced in 



V^ a G(x, P ,s) = [s+ P U(x)]'- a G(x,p, S ) 



In t space, 

Vl- a G(x,p,t) 



r(a) 



(13) 



— +pU{x) J (f _ r) i- a G(x,p,T)dT. (14) 



Thus, due to the long waiting times, the evolution of G(x,p, t) is non-Markovian 
and depends on the entire history. 
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In s space, the fractional Feynman-Kac equation reads 

d 2 

sG(x,p,s) - 5{x - x ) = Ka — ls+pUix)} 1 a G(x,p,s) - pU(x)G(x,p,s). 

(15) 

A few remarks should be made. 

(i) The integer Feynman-Kac equation. — As expected, for a — 1 our 
fractional equation (fT2")l reduces to the (integer) Feynman-Kac equation (JTJ. 

(ii) The fractional diffusion equation. — For p = 0, G(x,p — 0, t) = 
J Q G(x, A, t)dA reduces to G(x, t), the marginal PDF of finding the particle 
at x at time t regardless of the value of A. Correspondingly, Eq. (|10p reduces 
to the well-known Montroll- Weiss CTRW equation (for x = 0) [T2U55] : 



s 1 — cos(ka)ip(s) 
Eq. (|12p reduces to the fractional diffusion equation: 

_ G(:M)=i ^_^G(z,i), (16) 

where X'rl" is the Riemann-Liouville fractional derivative operator (2?j^L"G(:r. 
s 1 ~ a G(x, s) in Laplace f — > s space) [T2lf34] . 

(Hi) The scaling limit. — To derive our main result — the differential equa- 
tion (P|)— we used the scaling, or continuum, limit to CTRW [50115111521155] . 
In this limit, we take a — > and B a —> 0, but keep K a = a 2 /(2B a ) finite. 
Recently, trajectories of this process were shown to obey a certain class of 
stochastic Langevin equations 35 , 36 , 37] , hence giving these paths a mathe- 
matical meaning. 

(iv) How to solve the fractional Feynman-Kac equation. — To obtain the 
PDF of a functional A, the following recipe could be followed pQ: 

1. Solve Eq. (|15[) . the fractional Feynman-Kac equation in (x,p,s) space. 
Eq. (1151) is a second order, ordinary differential equation in x. 

2. Integrate the solution over all x to eliminate the dependence on the final 
position of the particle. 

3. Invert the solution (p, s) — > (A,t), to obtain G(A,t), the PDF of A at 
time t. 

We will later see (Section I2.2j) that the second step can be circumvented by 
using a backward equation. 

(v) A general functional. — When the functional is not necessarily posi- 
tive, the Laplace transform A — > p must be replaced by a Fourier transform. 
We show in the Appendix that in this case the fractional Feynman-Kac equa- 
tion looks like (|12l) . but with p replaced by —ip, 



^- t G(x,p,t)^K a ^V 1 t - a G(x,p,t)+i P U(x)G(x,p,t), (17) 

where G(x,p,t) is the Fourier transform A — > p of G(x,A,t) and T>]~ a — > 
[s — ipU(x)] 1 ^ a in Laplace s space. 
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(vi) Levy flights. — Consider CTRW with displacements A x distributed 
according to a symmetric PDF f{A x ) ~ |Z\j;|~( 1+AI \ with < \i < 2. For 
this distribution, the characteristic function is f(k) ~ 1 — C M |/c| M 12j. This 
process is known as a Levy flight, and as we show in the Appendix, the 
fractional Feynman-Kac equation for this case is (for A > 0) 

^G(x,p,t) = K atfi V%Vl- a G(x,p,t) -pU(x)G(x,p,t), (18) 

where K a ^ — C^/B a (units m^/sec"), and T>\~ a is the substantial fractional 
derivative operator defined above (Eqs. (fT3|) . (fl~4l) ) . V£ is the Riesz spatial 
fractional derivative operator defined in Fourier x — > k space as — > — {k^ 



2.2 A backward equation 



In many cases we are only interested in the distribution of the functional, A, 
regardless of the final position of the particle, x. Therefore, it turns out quite 
convenient (see Section [3]) to obtain an equation for G Xa (A, t) — the PDF of 
A at time t, given that the process has started at xq. 

According to the CTRW model, the particle, after its first jump at time 
t, is at either xo — a or xq + a. Alternatively, the particle does not move at 
all during the measurement time [0,t]. Hence, 



G X0 (A,t) 



j^(r) !±G X0+a [A 



TU(x ),t- t] 



+ -G X0 - a [A - tU{x ), t~r]jdr + W{t)8[A - tU(x )}. 

(19) 

Here, tU{xq) is the contribution to A from the pausing time at xq in the 
time interval [0, t\. The last term on the right hand side of Eq. (fTO)) describes 
a motionless particle, for which A(t) = tU(xo). We now Laplace transform 
Eq. P^|) with respect to A and t, using techniques similar to those used in 
the previous subsection. This leads to (for A > 0) 



G Xo {p,s) = -ip[s+pU(xo)] [G Xo+a (p,s) + G XQ - a (p,s)} 

l-$[s+pU(x Q )] 
s+pU(x ) 

Fourier transform xq —> ko of the last equation results in 

d 



G k {P,s) = Tp 



s+pU 



dkr 







s - 





cos(k Q a)G ko (p, s) 



S(k ). 
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As before, writing ip(s) ~ 1 — B a s a and cos(fc a) ~ 1 — a 2 k' 2 /2, we have 
d 



s + pU —i 



dkr 



G ko (p, s) + K a k G ko (p, s) 



6(k ), 



where we used the generalized diffusion coefficient K a = a 2 /(2B a ). Operat- 

r / \ i 1-a 

ing on both sides with 



sG ko (p,s) - 5(k ) = - K Q 



s + pU —i 



dkr 



ko 2 G ko (p,s) 







PU ( -i-Q)^ ) G ko (p,s). 



Inverting kr, — > xq and s — > t, we obtain the backward fractional Feynman- 
Kac equation: 



3 

g- t G X0 (p,t) = K a Vl- a -^G Xa {p,t) ~ P U(x )G Xo (p,t). (20) 

Here, T>]~ a equals in Laplace t — > s space [s +pU(xo)] 1 ~ a ■ The initial condi- 
tion is G Xo (A,t = 0) = 6(A), or G XQ (p,t = 0) = 1. In Eq. ([T2J the operators 
depend on x while in Eq. (|20p they depend on xq. Therefore, Eq. (TT2"j) is 
called the forward equation while Eq. ((20|) is called the backward equation. 
Notice that here, the fractional derivative operator appears to the left of the 
Laplacian d 2 /dxo 2 , in contrast to the forward equation (TT2"|) . 

In the general case when the functional is not necessarily positive and 
jumps are distributed according to a symmetric PDF f(A x ) ~ |Z\.j;|~( 1+A1 ), 
< /i < 2, the backward equation becomes (see the Appendix) 



G X0 ( P ,t) = K a ^V\~ a V^G X0 (p,t) + ipU(x Q )G X0 (p,t). (21) 



Here, p is the Fourier pair of A, T>\~ a — > [s — ipU(xo)] 1 ~ a in Laplace t — > s 
space, and V£ — > — |fco| M in Fourier x — > k space (see also comments (v) 
and (vi) at the end of section [2~T1 above) . 



3 Applications 

In this section, we describe a number of ways by which our equations can 
be solved to obtain the distribution, the moments, and other properties of 
functionals of interest. 
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3.1 Occupation time in half-space 

Define the occupation time of a particle in the positive half-space as T + = 
f*9[x(r)]dT (0(x) = 1 for x > and is zero otherwise) [THHH3]. To 
hnd the distribution of occupation times, we consider the backward equation 
(POj). transformed t -> s): 



s . r - ,„ ... , .^^^^(M * <0, 
^ a(s+p )i-«_|L Gxo( p )S) _ pG;co(P)S) ,,, .„. 



These are second order, ordinary differential equations in Xq. Solving the 
equations in each half-space separately, demanding that G Xo (p, s) is finite for 
|xo| -> oo, 

r ( n ,\- / C oexp(x W 2 /V^) + I x <0, 

Uxo(p ' s) -\c ie Kp[-x (s+pr/ 2 /viQ + ^ x >o. (Z6) 

For xq — > — oo, the particle is never at x > and thus G Xo (T + ,t) = 5(T + ) 
and G xo (p, s) = s , in accordance with Eq. (1231) . Similarly, for xo — > +oo, 
G Xo (T+,t) = 5(T+ -t) and G Xa (p,s) = (s+p) -1 . Demanding that G Xo (p,s) 
and its first derivative are continuous at xo = 0, we obtain a pair of equations 
forCo,C i: 

C + s- 1 = d + (s + p)- 1 ; C s Q/2 = -Ci(s+p) Q/2 , 

whose solution is 

^s+p)"/ 2 - 1 ps"/ 2 - 1 
— r— 75 , ; — ; — 77-777 I w 



[ s a/2 + ( s+p ) Q /2] ' * ( S+ p)[s«/2 + ( s+ p)«/2]' 

Assuming the process starts at xo = 0, Go(p, s) = G\ + (s +_p) _1 , or, after 
some simplifications: 

a/2-1 ,f , „W2-1 

<*<»'>- w^U^U - (24) 

Using [38 , the PDF of p + = for long times, is the (symmetric) Lamperti 

PDF: 

_ sin(7m/2) ( P+ )"/ 2 - 1 (l- P+ )"/ 2 - 1 

I[P+> 7T (p + )« + (l-p + ) a + 2(p + ) Q / 2 (l-p + )«/ 2 cos(7ra/2)' 

(25) 

This equation has been previously derived using different methods 16, 2"2"ll3"51 
[39] and was also shown to describe occupation times of on and off states in 
blinking quantum dots [TO1HT1H2"] . Naively, one expects the particle to spend 
about half the time at x > 0. In contrast, we learn from Eq. (|2"5"]) that the 
particle tends to spend most of the time at either x > or x < 0: f(p+) has 
two peaks at p+ = 1 and p + = (Fig. [IJ. This is exacerbated in the limit 
a — > 0, where the distribution converges to two delta functions at p+ = 1 and 
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at p+ = 0. For a = 1 (Brownian motion) we recover the well-known arcsine 
law of Levy PQUMISll^l ■ 

We note that the PDF (f2"5]) is a special case of the more general, two- 
parameter Lamperti PDF |16) : 

mp+ ) = «E&Mx (26) 

TT 

X (p+) a +K 2 (1- p+) a + 2K(p+) a / 2 (l - p+) a / 2 cos(7ra/2) ' 

where 1Z is the asymmetry parameter. In Eq. (|25l) . 7Z — 1 as a result of the 
symmetry of the walk. Consider, for example, the case when in Eq. ([2"2"| the 
diffusion coefficient is K£ for x < and K> for x > 0. Solving the equations 
as above, we obtain for f(p+) the two-parameter Lamperti distribution, Eq. 

HMD, with n = *Jk<Jk>. 

Kac proved in 1951 that for a — 1 (Markovian random-walk), the oc- 
cupation time distributions of both Brownian motion and Levy flights obey 
the same arcsine law [H]. It was therefore interesting to find out whether 
a similar statement holds for a < 1. We could not solve the Levy flights 
analog of Eq. (f2"2"j); therefore, we simulated trajectories whose PDF satisfies 
the fractional diffusion equation (Eq. {jX6j ) and its generalization to Levy 
flights (Eq. (Tig)) with p = 0). Simulations were performed using the subor- 
dination method described in [3"6ll37[|4"5] . The results are presented in Fig. [1] 
and demonstrate that indeed, for a < 1, the occupation time distribution is 
Lamperti's (j25|) for both /i = 2 and ijl < 2 (Levy flights). This result may be 
related to the recent finding that the first passage time distribution is also 
invariant to the value of /j, [46] , 



3.2 First passage time 

The time tf when a particle starting at xo = first hits x — b is called the 
first passage time and is a quantity subject to many studies in physics and 
other fields [37] • The distribution of first passage times for anomalous paths 
can be obtained from our fractional Feynman-Kac equation using an identity 
due to Kac [44] : 

Pi{t f >t} = Pr{ max x(t) <b}= lim G Xo (p,t), (27) 

0<T<f p— S-OO 

where the functional is Af — J Q U[x(t)]cIt, and 

. , | x < b, . „. 

^ = {l x>b. (28) 

This is true since G XQ (p, t) — J o °° e~ pA f G XQ (Af,t)dAf, and thus, if the parti- 
cle has never crossed x = b, we have Af = and e~ pAf = 1, while otherwise, 
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0.3 



L +0.2 



a=1 



O Q = 1,^ = 2 
□ a = 1, \x = 1.5 
a = l,M = l 
i3 a = 1, n = 0.5 

Lamperti, a = 1 (arcsinc) 

a = 0.5,// = 2 
V Q = 0.5, /i = 1.5 
< a = 0.5,^ = 1 
. q = 0.5,/i = 0.5 
— Lamperti, a = 0.5 




Fig. 1 The PDF of the occupation fractions in half- s pace T+/t . Trajectories of 
diffusing particles were generated using the methods of [36II37II45] with parameter 
values At = At = 10" 3 , At = 10" 2 (as defined in [35]), = 1, and x = 0. 
Simulations ended at t = 10 4 and include d 1 4 trajectories. Simulation results for 
a = 0.5, 1 and p, = 0.5, 1, 1.5, 2 (see Eq. (|21[) ~) are sho wn as symbols (see legend). 
Theoretical curves correspond to Lamperti's PDF, Eq. (|25[) (the arcsine distribution 
for a = 1), and are plotted as a solid line for a = 1 and as a dashed line for a = 0.5. 
It can be seen that the distribution of occupation fractions is determined by a but 
not by p.. 



Af > and for p — > oo, e pj4/ = 0. To find G Xo (p,t) we solve the following 
backward equation 



sG Xo (p,s) - I 



1-q 



.t < 6, 
(p, s) - pG^o (p, s) x > b. 



Solving these equations as in the previous subsection, demanding that G Xo (p, s) 
is finite for |a;o| — > oo and demanding continuity of G Xg (p,s) and its first 
derivative at xq = b, we obtain for xq = 



G (p,s) = - 



1 — e 



s a/2 _|_ ( s+p -j a /2 



To find the first passage time distribution we take the limit of infinite p, 

1 



1 / fc „a/2\ 

lim G (p, S ) = -(l-e-^ s ) 



(29) 



Defining r/ = (b 2 / Ka) l ' a , we invert s — > t: 

lim G (p,i) = Pr{t f > t} = 1 - / — Z Q 



/2 
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where I a /2(f) is the one-sided Levy distribution of order a/2, whose Laplace 

a/2 

transform is l a /2(s) = e~ s . The PDF of the first passage times, f(t), 
satisfies f(t) = §- t (Pr{t f < t}) = ^ (1 - Pr{</ > t}). Thus, 

/(*) = -lcc/2 (-) ■ (30) 

T f \ T fJ 

This result has been previously derived using different methods (e.g., Eq. 
(53) of [48]). The long times behavior of f(t) is obtained from the s — > 
limit: 



f(s) ~ 1 - 




Therefore, for long times 

fm-JFi^'^"'* (31) 

For a = 1, we reproduce the famous i~ 3 / 2 decay law of a one-dimensional 
random walk 47 . 



3.3 The maximal displacement 



The maximal displacement of a diffusing particle is a random variable whose 
study has been of recent interest (see, e.g., [4~9,50,51,52] and references 
therein). To obtain the distribution of this variable, we use the functional 
defined in the previous subsection (Eq. (1281 ). Let (t), and re- 

0<T<t 

call from Eq. (j2"T)) that Pr{x m < b} = lim G XQ (p, t). From the previous 

p— >oo 

subsection we have, for xq — (Eq. 



Pr{x m < b} = - ( 1 - e 7ft' 



Hence, the PDF of a; m is 

P(x m ,s) 
Inverting s t, we obtain 

P{%m 1 ^) 



t 



a\JK a 



s <x/2-l 

' 



>/2 



»/2 



1+2/a ' Q / 2 



2/a 



(32) 

This PDF has the same shape as the PDF of x up to a scale factor of 2 [3D] , 
and it is in agreement with the very recent result of |51) . derived using a 
renormalization group method. 
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3.4 The hitting probability 

The probability Ql{xo) of a particle starting at < xq < L to hit L before 
hitting is called the hitting (or exit) probability. The hitting probability 
has been investigated long time ago for Brownian particles |47j and more 
recently for some anomalous processes [S3]. For CTRW, it can be calculated 
using the following functional: 



U(x) 



< x < L, 
oo Otherwise. 



(33) 



With Eq. (|3"3"]) , A = L U[x(t)]cIt — as long as the particle did not leave the 
interval [0, L\ and is otherwise infinite. Therefore, G(x,p, t) — f °° e~ pA G(x, A, t)dA 
represents the probability of the particle to be at x at time t without ever 
leaving [0, L]. This is true for all p, since e~ pA is either or 1 regardless of p. 
At the boundaries, G(x = Q,p, t) = G(x — L,p, t) = 0. At (0, L), the forward 
fractional Feynman-Kac equation (Eq. (|15|1 ) reads, in s space, 



d 2 

sG(x,s) - S(x - xo) = Ka-s 1 a —^rG(x,s). 

ox 



(34) 



Note that Eq. ([34|l does not depend on p and is equivalent to the fractional 
diffusion equation, Eq. (|16p. with absorbing boundary conditions. The solu- 
tion of Eq. for x ^ xq is 



G(x, s) 






X 











(L - x) 



X < Xo, 
X > Xq- 



Matching the solution at x = xo and demanding -^G(x = Xq, s) — ^G(a 
x^, s) = — K ] ji _ a (from Eq. (J2U)), we have, for x > xq, 



G{x,s) 



1 



sinh 



sinh 



1/2 



sinh 



s «/2 



-.(L-x) 



>x - (35) 



The flux of particles that have never before left [0,L] and that are leaving 
[0,i] at time t through the right boundary is [53] 

J(L,t) = -K a V)^ a t ^-G{x = L,t), 

where £>rl" 1S the Riemann-Liouville fractional derivative, equal to s 1_Q in 
Laplace t — > s space (see Eq. ([TS])'). The hitting probability is the sum over 
all times of the flux through L [37] : 





Ql(xo)= / J{L,t)dt 
Jo 



K a s L - a —G(x = L,s) 
ox 



s=0 



14 



Using Eq. (|35|) . we have 

Ql(x q ) = j ; . (36) 

The hitting probability for anomalous diffusion, a < 1, is the same as in the 
Brownian case |47j . This is expected, since the hitting probability should not 
depend on the waiting time PDF ip(r). 

Note that a backward equation for Ql(xq) can be obtained by the much 
simpler argument that for unbiased CTRW on a lattice, Ql(xo) = [Ql(%o + a ) + Ql(xo 
In the continuum limit, a — > 0, this gives d ^s" * 1 = 0. With the boundary 

conditions Ql{xq = 0) = and Ql(xq = L) = 1, Eq. (fMf immediately 
follows (see |47j for a binomial random walk) . 



3.5 The time in an interval 

Consider the time-in-interval functional Ti = J* U[x(t)]cIt, where 



U(x) = 



1 \x\ < 6, 
|x| > b. 



(37) 



Namely, Tj is the total residence time of the particle in the interval [—6,6]. 
Denote by G Xo (Ti,t) the PDF of T$ at time t when the process starts at xq, 
and denote by G Xo (p,s) the Laplace transform 7$ — > p, t — > s of G Xo (Ti,t). 
G Xo (p,s) satishes the backward fractional Feynman-Kac equation: 



sG Xo (p,s) 



1 



K a {s+p) 1 a -§^G Xa {p,s) -pG Xo (p,s) \x \ < 6, 



K^-v^GxJp, s) 



dx 2~X 0K i>1»J l^0| > 

We solve this equation demanding that the solution is finite for |xq| — > oo 



(38) 



G Xo (p, s) 



_ fc lC osh[xo(s+p) a / 2 /VTC\ +^ M<6, 



C exp[-|x |W 2 /V7^ 



+p '"- u '^"' (39) 
x | > b. 



Demanding continuity of G Xo (p, s) and its first derivative at xq = b we solve 
for C\ and then obtain for xq = 



Ga(p, s) = 



P 



cosh 



\/kZ 



,/2 



sinh 



(s+P) 



</2 



»(»+P){ 



cosh 



^f^sinh 



(40) 



In principle, the PDF Go (Ti , t) can be obtained from (j40|) by inverse Laplace 
transforming p — > and s — > t. However, we could invert Eq. f)40[) only for 
a -> 0: 



GoCTi,*)^ = (1 - e^^)^ - 1) + e- 6 /^(Ti). 



(41) 



This can be intuitively explained as follows. For a — > 0, the PDF of a; becomes 
time-independent and approaches G(x,t) ss exp(— \x\ / t/Kq) / (2t/Kq) (Eq. 
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(Al) in [30 ). With probability j b _ b G{x,t)dx = 1 - e ~ b /^, the particle 

never leaves the region [—6, b] and thus Ti = with probability e~ h l*^° , the 
particle is almost never at [—6, 6] and thus Ti = 0. 

The first two moments of Ti can be obtained from Eq. (|40[) by 



(Ti)(s) = -^G ( P ,s) 



(T?) (.) 



dp 



2 G o(P,s) 



p=0 



Calculating the derivatives, substituting p 
the long times limit, 



0, and inverting, we obtain, in 



(Ti 



l-a/2_ 



/K^r(2-a/2)' 



26(1 - a) 



/E^r(3 - a/2) 



;-a b 2 (3a-l) 
K a r(3 - a) ' 



(42) 



We verified that Eq. (|42p agrees with simulations (Fig. [5]) . The average time 
at [—6, 6] scales as i 1- "/ 2 since this is the product of the average number of 
returns to the interval [—6, b] (~ W 2 ) and the average time spent at [—b, b] 
on each visit (~ i 1_Q ; see Eq. (61) in [IB]). We also see that for a < 1, the 
PDF of T 4 cannot have a scaling form since <T 2 ) - ^ ^2 _ f2 -a 

For a = 1, (Ti) ~ i 1 / 2 and <T 2 ) - i. 



3.6 Survival in a medium with an absorbing interval 

A problem related to that of the previous subsection is a medium in which a 
diffusing particle is absorbed at rate R whenever it is in the interval [—b, b]. 
The survival probability of the particle, S 1 , is related to Ti, the total time at 
[-6, b], through S = exp{-RTi). Thus, if G Xo (T u t) is the PDF of T, at time 
t, then the Ti — > R Laplace transform G XQ (R,t) = J" 00 e~ RT *G Xo (Ti,t)dTi 
equals (S), the survival probability averaged over all trajectories [21]. From 
Eq. (j4T)]) of the previous subsection we immediately obtain (in Laplace t — > s 
space and for xq = 0) 



(S) = G (R,s) 



R+si cosh I {s+ S- /2 b 



{s+R) c 



■ sinh 



(s+nr /2 , 



3(S + fl) { 



cosh 



(s+fl)"/ 2 z. 



(s+fl) c 



■ sinh 



(43) 

where here R is a parameter (the absorption rate) and thus the equation 
needs to be inverted only with respect to s. We could invert (H31 for a few 
limiting cases. 

(i) t — > oo. The long time behavior is obtained by taking the s — > limit 
and inverting: 



(S) 



r(l - a/2) sinh 



bR a ' 2 



(Rt) 



-a/2 



0[(Rt)- a ] 



(44) 
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Fig. 2 The first two moments of the time-in-interval functional Ti. Top panel- 
first moment; bottom pan el- second moment. Simulations (circles) were performed 
using the method of [36II37] with parameters as i n Fig. [1] The theoretical curves 
(line s) correspond to the leading terms of Eq. Q42p . The width of the interval (Eq. 
p7)) is b = 0.1. 

Thus, the survival probability of the particle in the absorbing domain decays 
as t~ a l 2 . We verified Eq. (l44l) using simulations (Fig. [3]). 
(ii) a — > 0. Inverting Eq. (|43p yields 

(S) a ^ Q = (1 - e -i/VKJ) e -» + e -6/VKo. (45) 
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Fig. 3 Survival probability in an absorbing medium. We simulated anomalous 
diffusion trajectories as in Fig. fTI [361157] . with total time t = 10 3 . We plot (e~ RTl ), 
where Ti is the total time at [—6,6], R = 10, and 6=1 (Eq. (|37[0 . Simulation 
results (symbols) agree with theory (lines), Eq. (|44|) . for Rt 3> 1. 



This can be explained as in the previous subsection. For a — ¥ 0, the PDF of a; 
approaches G(x,t) « exp(-\x\/^/K^)/(2^/K^). With probability (l - er b l^ 
the particle never leaves the region [—6,6]. Thus, its probability of survival 
is just e~ Rt . With probability e~ h l^~° \ the particle is almost never in the 
absorbing zone, and it survives with probability 1. 

(Hi) Other limiting cases. It can be shown that for 6 — > or R — > 0, 
(S) = 1; for R^oo, (S) = 0; and for 6 -> oo, (S) = e~ m . 



3.7 The area under the random walk curve 

The functional A x = J * x(t)<1t (U(x) = x) represents the total area under 
the random walk curve x(t) [51125]. and it is also related to the phase accumu- 
lated by spins in an NMR experiment [35J. In this subsection we obtain the 
first two moments of this functional, and for a couple of special cases, also 
its PDF. Since A x is not necessarily positive, we use the generalized forward 
equation (Eq. (TT7]) Laplace transformed t —> s), 

d 2 

sG(x,p 7 s) - S(x) = K a -^--^(s - ipx) 1 a G(x,p,s) + ipxG(x,p,s). (46) 

Here, G(x,p,s) is the Fourier-Laplace transform of G(x,A x ,t) and we as- 
sumed xq = 0. Since the walk is unbiased, (A x ) = 0. To find the second 
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moment of A, T , wc use 



<4> (*) = 



52 



2 G(x,p,t) 



dx. 



p=0 



Integrating Eq. (|4"6")l over all x, taking the derivatives with respect to p and 
substituting p = 0, we obtain 



,.<4>( S ) = 2^ X )( S ), 



(47) 



which is in fact obvious since (A x ) = ^J* 'x(T)<iT^ = x, and thus 



± / A2 
dt 



{A 2 ,} = 2(xA x ). Hence, the problem of finding \A X \ reduces to that of 



dx. This 



p=0 



finding (xA x ), for which we have (xA x ) — —ix -^G(x,p,t) 
leads to 

S (xA x )( s ) = (x 2 }(s). (48) 

Similarly, 



:(x 2 ) (s) = 2K aS - a . 



(49) 



Combining Eqs. 071, 05J|, and (gU), we find (s) = AK aS -^ +a \ or, in 

f space, 



(A 2 x )(t) 



4K« _ t 2 +a 



r(3 + a) 



(50) 



Higher moments of A x can be similarly calculated (see next subsection) . The 
distribution of A x can be obtained for a few limiting cases. For a = 1, A x is 
normally distributed (Eq. (61) in [S]): 



G(A«,t) £ 



471-^1 £ 3 



exp 



(51) 



For a -> 0, the PDF of x is G(^,i) w exp(-|a;|/v A Ko)/(2v A Ko) ([30] and 
Section [3.5[) and is independent of t. In other words, the particle is found at 
x(t) for most of the time interval [0,t]. Hence, A x (t) ~ to(t) and 



1 



exp 



I A, 



(52) 



To confirm Eqs. (|51l) and (I5"2"j) . we plot in Fig.0]the PDF of A x for various 
values of a as obtained from simulation of diffusion trajectories. It can also 
be seen from Fig. [4] that the PDF of A x obeys a scaling relation, as we show 
in the next subsection. 
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Fig. 4 The area under the random walk curve A x . We simulated diffusion trajec- 
tories as in Fig. □ [361137] (10 5 trajectories) and calculated the PDF G{A x ,t). To 
illustrate the scaling, we plot t 1+a/2 G(A x , i) vs. A x /t 1+a/2 (K a = 1), collapsing 
all curves with the same a but different times: t — 1 (circ les), t — 10 (squares), 
and t = 100 (diamonds). The ory for a = 1 is from Eq. (|51|l (Gaussian, solid line); 
theory for a — > is from Eq. (|52l) (exponential, dashed line). 



3.8 The moments of the functionals U(x) = x k 

In the previous subsection we derived the first two moments of the U{x) = x 
functional; but in fact, all moments of all functionals A x h = J Q x k {r)dT, 
k = 1,2,3, ... can be obtained, leading to a scaling form of their PDF. As 
explained above, the functionals with k = 1,2 arise in the context of NMR 
and are therefore particularly interesting. 

We assume xq — and consider the forward equation (fl~5|) for even fe's: 



sG(x,p,s) —S(x) = K a -g-^ (s 



px J G{x,p,s) 



px k G(x,p,s). (53) 



Here, G(x,p,s) is the double Laplace transform of G(x,A x k,t) since for 
even k's A x k is always positive. We are interested in the moments (A™ k }, 
n = 0, 1, 2, ...; however, to find these, we must first obtain the more general 
moments (A n x m ), n,m = 0, 1, 2, .... Operating on each term of Eq. (|53"|) with 
(— l) n J^r , substituting p = 0, multiplying each term by x m , and integrating 
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over all x, Eq. ([53]) becomes 

a (A n x m ) (s) = 8 n , 5 m , + H n _ x n (A^x^) (s) + 

(-iy s 1 - 01 ^ x 

x (A^x 171 ^- 2 ) (s), (54) 

where 5ij is Kronecker's delta function — 5ij equals 1 for i = j and equals 
zero otherwise; and Hi is the discrete Heaviside function — Hi equals 1 for 
i > and equals zero otherwise. It can be proved that Eq. ([54"]) remains true 
also for odd fc's, when A x k can be either positive or negative. Eq. ([54]) is 
satisfied by the following choice of (A n x rn ): 

(A n x m ) (s) = c, lhm (k)K^ s -(i+n+^ a ) i (55) 



H m ^ 2 K a m(m - 1) 



E 

.7=0 



3-1 
.1=0 



a -I) 



for all n and even m when k is even and for even (n + m) when k is odd. 
In all other cases (A n x m ) = due to symmetry. The c„. m 's are fc-dependent 
dimensionless constants that satisfy the following recursion equation: 

(k) = SnftSmfi + Hn-1 (*)+ 

(-l) J c n _ :)> „ +:)fc _2(fc), (56) 



H m - 2 m(m - 1) 



7=0 



J=0 



a-/) 



with initial conditions co,o(&) — 1 and Co,i(A;) = 0. The moments of A x k are 
therefore given in t space by 



(A n xk ) (t) = c „, (fc) 



r i 



nak 
2 



(57) 



For example, for jfe = 1, = = 0, = 4 J ftT Q i 2 + a / / r(3 + a) (Eq. 

([50])). and (A*) = 48(a 2 + 7a + 12)^i 4 + 2Q / J T(5 + 2a); while for fc = 2, 
(A^) = 2K a t 1+a /r{2 + a) and (A 2 x2 ) = (48 + 8a)i^ 2 + 2a /r(3 + 2a). 
Eq. ([57]) suggests that the PDF of A,.* obeys the scaling relation 

G(Axk ' t] = K *A+<*i* 9a ' k (jf*/v^/>) ' (58) 

where g a ,k(x) is a dimensionless scaling function. To verify the scaling form 
of Eq. ([58]) . we plot in Fig. @] simulation results for the PDF of A x (k = 1) 
for a ~ and a = 1 (for which G(A x ,t) is known — Eqs. ([5Tj) and (1521 in 
the previous subsection), and for an intermediate value, a = 0.5. In all cases 
the simulated PDF satisfies the scaling form (|58|). 
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4 Summary and discussion 

Functionate of the path of a Brownian particle have been investigated in nu- 
merous studies since the development of the Feynman-Kac equation in 1949. 
However, an analog equation for functionate of non-Brownian particles has 
been missing. Here, we developed such an equation based on the CTRW 
model with broadly distributed waiting times. We derived forward and back- 
ward equations (Eqs. (fT^j) and ([20)) ) and generalizations to Levy flights (Eqs. 
(UHJ) and (|2"T)) ). Using the backward equation, we derived the PDFs of the 
occupation time in half-space, the first passage time, and the maximal dis- 
placement, and calculated the average survival probability in an absorbing 
medium. Using the forward equation, we calculated the hitting probability 
and all the moments of U(x) = x k functionate. 

The fractional Feynman-Kac equation (|12[) can be obtained from the in- 
teger equation ((TJ) by insertion of a substantial fractional derivative operator 
[25] , In that sense, our work is a natural generalization of that of Kac's. The 
distributions we obtained for specific functionate are also the expected ex- 
tensions of their Brownian counterparts: the arcsine law for the occupation 
time in half-space [Tll43j was replaced by Lamperti's PDF (Eq. ((25)) ) 22J, 
and the famous i~ 3 / 2 decay of the one-dimensional first passage time PDF 
[47] became t~( 1+Q / 2 ) (Eq. ([5T10 . Thus, our analysis supports the notion that 
CTRW and the emerging fractional paths 36, 37 are elegant generalizations 
of ordinary Brownian motion. Nevertheless, other non-Brownian processes 
are also important. For example, it would be interesting to find an equa- 
tion for the PDF of anomalous functionate when the underlying process is 
fractional Brownian motion |14j . 

Our fractional Feynman-Kac equation (|T2)) has the form of a fractional 
Schrodinger equation in imaginary time. Real time, fractional Schrodinger 
equations for the wave function have also been recently proposed 55,56,57, 
58,59 . However, these are very different from our fractional Feynman-Kac 
equation. In 55,56,57 , the Laplacian was replaced with a fractional spatial 
derivative which would correspond to a Markovian CTRW with heavy tailed 
distribution of jump lengths (Levy flights; see also the Appendix below). 
The approach in [58,59^ is based on a temporal fractional Riemann-Liouville 
derivative — however not substantial — which leads to non-Hermitian evo- 
lution and hence non-normalizable quantum mechanics. It is unclear yet 
whether all these fractional Schrodinger equations actually describe any phys- 
ical phenomenon (see |60j for discussion) . In principle, a fractional Schrodinger 
equation can also be written using the substantial fractional derivative we 
used here. If there is a physical process behind such a quantum mechanical 
analog of our equation remains at this stage unclear. 

In this paper we considered only the case of a free particle. In [29] . we 
reported a fractional Feynman-Kac equation for a particle under the influence 
of a binding force, where anomalous diffusion can lead to weak ergodicity 
breaking [rJTl[rj2llB3] . The derivation of an equation for the distribution of 
general functionate and the treatment of specific functionate for bounded 
particles will be published elsewhere. 
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Appendix: Generalization to arbitrary functionals and Levy flights 



Here we generalize our forward and backward fractional Fcynman-Kac equa- 
tions f dl^j) and (PH)) . respectively) to the case when the functional is not 
necessarily positive and to the case when the CTRW jump length distribu- 
tion is arbitrary, and in particular, heavy tailed. 

In our generalized CTRW model, the particle moves, after waiting at x, 
to x + A x , where A x is distributed according to f(A x ). The PDF f(A x ) 
must be symmetric: f(A x ) = f(—A x ) but can be otherwise arbitrary. Let us 
rederive the forward equation for this model. We replace Eq. Q with 

(•t fOO 

Q n+1 (x,A,t)= 4>{t) f(A x )Q n [x- A x ,A-tU(x- A x ) : t-T]dA x dr. 

JO J -oo 

Since A can be negative, we Fourier transform the last equation A — > p 

ft fOO 

Q n+1 (x,p,t) = / ^(r) / f(A x )e lpTU( - x - A ^Q n (x-A x ,p,t-r)dA x dT. 

JO J -oo 

Laplace transforming t — > s and Fourier transforming x — > k we have 



sn+i 



(k,p, s) 



ikx 



f(A x )iij)[s-ipU(x-A x )]Q n (x-A x ,p, s)dA x dx. 



Changing variables: x' = x — A x , 



Q n+ i(k,p,s) 



i k _A . 



f(A x )dA a 



/(fc)V 



s-i P U\-i- 



e tkx 'i)\s - ipU(x')]Q n (x',p, s)dx' 

Qn(k,p, s). 



Summing over all n and using the initial condition Qo(k,p, s) 



J2Qn(k,P,s) = ll-f(k)i> 
n v 



s — ipU 



-4) 



o ikx 



Note that this agrees with Eq. ([5} since for nearest neighbor hopping f(k) = 
f^oo e lfczix [§£(Ar - a) + \5{A X + a)] dA x = cos(fca). Next, we observe that 
Eq. © of Section I!Q1 remains the same even under the general conditions. 
Calculating the transformed G(k,p,s) as above, and using the result of the 
last equation, we obtain the formal solution 



G(k,p,s) 



.i-V^-'WHm)] 



s-ipU(-i^) 



s-ipUl-i- 



ikxQ 
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We now assume that f{A x ) has a hnite second moment and thus its char- 
acteristic function can be written, for small k, as f(k) ~ 1 — <r 2 k 2 /2. This 
characteristic function is identical to that of nearest neighbor hopping (with 
a = a); we can thus proceed as in Section [2. II to obtain 

^-G(x,p,t)=K a -^V 1 t - a G(x,p,t)+i P U(x)G(x,p,t), (59) 
at ox z 

where here T)\~ a — > [s — ipU (a;)] 1- " in Laplace s space and K a = a 2 / (2B a ). 

Consider now the case of Levy flights — f(A x ) ~ |zA a ;| _ ( 1+A ') (for large 
A x ) with < p < 2, and thus jump lengths have a diverging second moment. 
The characteristic function is f(k) ~ 1 — C A ,|fc| M , and the fractional Feynman- 
Kac equation becomes 

^ t G(x,p,t)=K a ^V 1 t - a G(x,p,t) + ipU(x)G(x,p,t), (60) 

where K a ^ = C^j B a and is the Riesz spatial fractional derivative oper- 
ator: — > — |fc| M in Fourier k space. 

Repeating the calculations of Section 12.21 for a non-necessarily-positive 
functional and for Levy flights, it can be shown that the generalized backward 
equation is: 

^G X0 (p,t) = K a ^Vl- a V% G Xo (p,t) + ipU(x )G Xo (p,t). (61) 

Here, T>\~ a — > [s — ipU(xo)} 1 ^ 01 in Laplace s space and V£ — > — l^'ol^ in 
Fourier kg space. 
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